using namespace std;

//allow use of global variables from table
extern vector<double> tab_mass_axis;
extern vector<double> tab_low_mass_axis;
//extern vector<double> tab_low_nuv_slope;
//extern vector<double> tab_low_fuv_slope;
//extern vector<double> tab_low_bol_slope;

//added below vectors (RdS)
//extern vector<double> tab_low_bol;
//extern vector<double> tab_low_fuv;
//extern vector<double> tab_low_nuv;
//extern double** low_matrix;

extern long ntables;
extern double (*imf)(gsl_rng *);
extern vector<double> low_fluxes;
extern double kludge_index;

/* procedure that populate the cluster according to a given mass function */

double star_draw(double cluster_mass, double evolve_limitmass,
		       vector<double> *star_masses, vector<double> *star_indexes,
			gsl_rng * rand, double *kcount, double *biggest_star, 
			double ** low_matrix){

   
  //populate that cluster until it is full, randomly
  double total_mass=0;
  double new_mass=0;
  double lowmass_index;  //not follow low mass stars
  long floor_indx=0;
  long ceil_indx=0;
  double mod_indx=0;
//  double kludge_index=20.98;////get_tab_indexes(tab_low_mass_axis,0.151);
//  double sub_mass=0;//this is a holder for value of the last cluster if there is
			//overfill

  *kcount=(double) 0;
  *biggest_star=0;
  //float temp_mass_ind=0;
//  while (total_mass < cluster_mass){
    //pick randomly from the imf
    new_mass=(*imf)(rand);

    total_mass+=new_mass;
    
    //now update arrays, if adding the star gets you closer to the drawn mass
    //then do so... if not then don't
//   if (fabs(total_mass-cluster_mass)  < fabs((total_mass-new_mass)-cluster_mass)){  
							//check such that if it 
							//is closer to include the star then
							//do so
      if (new_mass > evolve_limitmass){//non-low mass star case
	//store the mass
	star_masses->push_back(new_mass);
      	if (new_mass > *biggest_star){//check if biggest star... if so then save
	*biggest_star=new_mass;
      }
      } else {//low-mass star case
	

	//Get the indexes
	//temporary fix for the fact that the table doesn't go below 0.151 M
	//If so, just redefine the mass to be 0.151 (do not change cluster mass)
	if(new_mass>=0.156) {
	  lowmass_index=get_tab_indexes_low(tab_low_mass_axis,new_mass); 
	}  else {
	  lowmass_index=kludge_index;      //get_tab_indexes(tab_low_mass_axis,0.151);
					//our tables don't go below 0.151 solar masses
					//this number has to be updated for a new table
	}
	//fill in the luminosity and forget about them
	//first find index for this star for low-mass tables
	floor_indx=(long) floor(lowmass_index);
	ceil_indx=(long) ceil(lowmass_index);
	mod_indx=lowmass_index-floor_indx;

	for (long ntab=0;ntab<ntables;ntab++){//calculate all desired fluxes and add
	low_fluxes[ntab]+=pow(10,(low_matrix[ntab][ceil_indx]-
		low_matrix[ntab][floor_indx])*mod_indx+low_matrix[ntab][floor_indx]);
	  }
	new_mass=1e5; //this is to make the star have a high deathmass so its
			//properties are not evaluated
	       star_masses->push_back(new_mass);
      }//end low mass star case
      
    *kcount=*kcount+1;//this counter keeps track of number of all stars
//    } else sub_mass=new_mass;//if the cluster is overfilled then subtract out last iteration
    
//  }//end while to fill up cluster

  //get the real total mass (if overfilled)
  //total_mass-=sub_mass; //submass=0 unless there is overfilling

  return total_mass;
  
}







